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ABSTRACT 


An  implicit  finite  difference  version  of  the  15°  parabolic  equation 
first  developed  by  Claerbout  was  used  to  model  acoustic  wave  propagation 
in  shallow  water.  The  algorithm  uses  a  variable  grid  spacing  in  the 
depth  as  well  as  range  direction,  resulting  in  rapid  execution.  The 
water-sediment  interface  was  simulated  by  an  Epstein  layer.  An  attempt 
to  model  propagation  loss  data  on  the  continental  slope  and  shelf  near 
Nova  Scotia  was  unsuccessful  because  of  a  lack  of  adequate  environmental 
data,  as  well  as  deficiencies  in  the  modeling  method.  However,  from 
modeling  it  was  found  that  the  sediment  properties  controlled 
propagation  loss  near  the  cutoff  frequency. 
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I.  INTRODUCTION 

Models  for  the  propagation  of  underwater  sound  in  shallow  water  are 
complicated  by  the  interaction  with  the  ocean  bottom.  When  the 
propagation  is  range— dependent ,  the  models  become  even  more  complex. 
One  method  which  has  been  investigated  as  a  tool  for  analyzing  shallow 
water  problems  is  the  parabolic  equation  technique  [Lee  and  Gilbert, 
1982].  Since  the  parabolic  equation  is  itself  an  approximation  to  the 
wave  equation,  however,  special  techniques  must  be  used  to  model 
acoustic  interfaces.  This  thesis  is  devoted  to  the  development  of  a 
parabolic  equation  method  for  modeling  propagation  loss  in  a  range- 
dependent  shallow  water  environment.  The  method  is  tested  by  attempting 
to  model  actual  data  recorded  in  two  separate  shallow  water 
environments . 

Explosive  signals  generated  on  the  continental  shelf  and  slope  near 
Nova  Scotia  [Brocher  et  al . ,  1981]  were  analyzed  to  study  acoustic 
propagation  in  a  range-dependent  shallow  water  environment.  The  data 
were  recorded  by  ocean  bottom  sensors  at  two  water  depths.  The  shallow 
sensor,  at  a  depth  of  67  m,  was  situated  in  a  region  of  bottom-limited 
sound  propagation.  The  other  sensor  was  on  the  continental  slope  in 
1301  m  of  water.  Seismic  propagation  loss  measurements  were  obtained 
from  the  recorded  signals  [Brocher  et  al.,  1982].  With  geoacoustic 
models  of  the  water  column  and  sediment  structure  obtained  from  various 
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sources,  the  parabolic  equation  was  used  to  model  the  observed 
propagation  loss* 

The  parabolic  equation  is  a  useful  means  of  calculating  long  range 
low  frequency  acoustic  propagation  when  the  raypath  turning  points  occur 
within  the  water  column  [Hanna,  1976;  Hanna  and  Rost,  1981]#  This  type 
of  propagation  usually  occurs  in  a  deep  ocean  with  a  well  defined  sound 
channel*  In  shallower  areas,  where  bottom  interacting  paths  are 
present,  the  parabolic  equation  method  has  been  less  successful  [Lee  and 
Gilbert,  1982].  As  a  result,  there  is  a  lack  of  literature  describing 
parabolic  equation  modeling  of  actual  shallow  water  propagation  loss 
data*  The  difficulty  stems  from  the  discontinuous  change  in  velocity 
and  density  at  the  bottom  of  the  ocean*  Frisk  et  al  •  [1981]  derived  a 
parabolic  equation  technique  that  accurately  modeled  the  water-sediment 
interface  for  single  bottom  bounce  paths  in  the  North  Atlantic,  but 
this  technique  is  not  applicable  to  mul t i— bounce  paths.  This  thesis 
presents  a  more  general  method  of  treating  acoustic  interfaces  for 
calculations  based  on  the  parabolic  equation. 

Normal  mode  theory  gives  an  exact  solution  to  the  wave  equation  in 
all  horizontally  stratified  media.  However,  modifying  the  normal  mode 
technique  to  take  into  account  range-dependent  variation  in  acoustic 
media  has  met  with  limited  success  [Graves  et  al.,  1975].  In  contrast 
with  normal  modes,  the  parabolic  equation  is  derived  from  the  wave 
equation  without  recourse  to  separation  of  variables.  Thus,  the 
parabolic  equation  technique  is  not  restricted  to  the  range-independent 
environment.  For  the  same  reason,  mode  coupling,  a  phenomenon 
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associated  with  sloping  bottoms  and  other  two  dimensional  media,  is 
handled  by  the  parabolic  equation  [Jensen  and  Kuperman,  1980]*  In  spite 
of  its  difficulty  with  interfaces,  the  parabolic  equation's  capability 
to  solve  range-dependent  problems  and  its  correct  modeling  of  mode 
coupling  make  it  a  natural  candidate  for  modeling  shallow  water  acoustic 
propagation, 

McDaniel  [1975b]  and  Lee  et  al,  [1981]  devised  implicit  finite 
difference  algorithms  to  solve  the  parabolic  equation.  Another  method 
introduced  by  Tappert  and  Hardin  combines  the  use  of  the  Fourier 
transform  with  a  finite-difference  operation  [Tappert,  1977],  Of  these 
two  techniques,  the  implicit  finite  difference  method  was  selected  for 
the  modeling  reported  here  because  it  can  be  used  with  a  variety  of 
boundary  conditions  [DiNapoli  and  Deavenport,  1979]*  Another  advantage 
of  implicit  finite  differences  is  faster  execution,  McDaniel  [  1975b] 
found  this  method  to  be  significantly  faster  than  that  of  Tappert  and 
Hardin.  Here,  we  use  an  implicit  finite  difference  algorithm  to  model 
acoustic  wave  propagation.  A  variable  grid  spacing  is  used  that  permits 
propagation  loss  models  to  be  computed  in  less  time. 
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II.  DEVELOPMENT  OF  NUMERICAL  MODELING  TECHNIQUES 


Derivation  of  the  parabolic  approximation  to  the  wave  equation 


The  wave  equation  using  the  Laplacian  operator V 
and  the  field  variable  being  pressure  p  is 


2 


with  velocity  c, 


v2P  -  Ip  -  o  .  [i] 

2 

c 

Inserting  a  continuous  sinusoidal  time  dependence  and  changing  to 
cylindrical  coordinates  with  no  azimuthal  dependence  results  in  the 
Helmholtz  or  reduced  wave  equation, 


rr 


1 

-P, 


zz 


k2p  =  0 


[2] 


r 

where  k  =  k(z,r),  c  =  c(z,r),  k  =  o)/c,  u>  *  2irf.  The  range,  r,  is  in  the 
principal  direction  of  propagation  and  z  is  depth,  positive  downward. 
Since  equation  [2]  is  an  elliptic  partial  differential  equation  with 
second  derivatives  in  both  r  and  z,  it  is  difficult  to  solve 
numerically.  Some  approximations  will  be  made  in  order  to  reduce 
equation  [2]  to  a  simpler  form.  In  addition,  two  changes  of  variables 
will  serve  to  decrease  the  spatial  variation  of  the  computed  wavefield, 
resulting  in  less  numerical  error. 
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It  is  usual  to  normalize  for  cylindrical  spreading  by  introducing 


P  = 


[3] 


which  gives 


v 


+  v 


[4] 


rr  zz  7 

4r 

Next,  a  far-field  approximation  can  be  made  by  neglecting  the  term 

9  9  2 

-l/(4r  )  since  v/r  will  be  much  smaller  than  k  v  for  ranges  on  the 

order  of  a  few  wavelengths  or  greater.  Using  a  "square-root  operator 

and  binomial  expansion  [Claerbout,  1970b]  or  splitting  matrix  [McDaniel, 

1975b],  equation  [4]  may  be  reduced  to 


v  =  i(k  + _ )v  .  [5] 

r  2 
2k  dz 

The  solution  to  this  parabolic  equation  represents  an  outward  traveling 
wave  field.  The  loss  of  the  wavefield  second  derivative  with  respect  to 
r  forces  a  restriction  on  field  gradients  that  can  be  treated  using 
this  approximation.  The  equation  of  the  inward  traveling  wave  field  has 
been  decoupled  from  equation  [5].  In  the  splitting  matrix  technique,  a 
matrix  equation  arises  with  terms  that  cause  coupling  between  the 
outward  and  inward  traveling  waves  [Corones,  1975].  These  terms  are 
then  explicitly  deleted  to  give  equation  [5] •  Physically,  coupling 
manifests  itself  as  back-reflection. 
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A  constant  "average"  wavenumber  =  io/cq  is  chosen  and  the 
following  change  of  variable  is  used 

v  ■  u  exp  [  ikpr  ]  [  6  ] 

The  new  field  variable  u  changes  slowly  with  range.  This  result  allows 
larger  step  sizes  in  the  numerical  implementation  of  the  parabolic 
equation.  The  effect  of  the  substitution  is  like  moving  the  coordinate 
system  at  the  constant  velocity  Cq  =  to / k^  [Volk,  1975].  In  practice,  kg 
is  selected  so  that  Cq  is  near  the  average  velocity  of  the  acoustic 
medium.  Inserting  equation  [6]  into  equation  [5]  gives  the  final 
parabolic  equation, 

u  =  (  i(k  -  kn)  +  ^ _ )u  ,  [7] 

r  0  2 

2k  dz 

which  is  the  basis  for  the  computer  model  presented  in  this  thesis. 

The  following  comment a  pertain  to  both  equation  [7]  and  the  similar 
Tappert  and  Hardin  [Tappert,  1977]  equation.  The  parabolic  equation  [7] 
gives  an  accurate  solution  for  the  propagation  of  a  single  normal  mode 
with  wavenumber  k^  if  the  waveguide  does  not  induce  coupling  into  other 
modes  [Fitzgerald,  1975] .  For  modes  of  propagation  with  values  of 
wavenumber  k  departing  from  k^  (or  phase  velocity  c  differing  from  Cq) 
the  wavefield  contains  increasing  errors  in  the  phase  velocities  and 
group  velocities  of  propagating  modes  [Claerbout,  1970b].  Equation  [7] 
is  thus  known  as  the  15°  or  narrow  bandwidth  approximation  since  it 
accurately  propagates  a  cone  of  rays  with  a  spread  of  about  15°  or  less. 
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Dropping  the  second  derivative  urr  leads  to  greater  errors  for 
nonhorizontal  raypaths  [Volk,  1975]  and  also  requires  that 
discontinuities  in  acoustic  properties,  as  at  the  ocean-sediment 
interface,  must  be  smoothed  in  order  to  apply  the  method  [Tappert, 
1977]  . 

Fitzgerald  [1975]  compared  parabolic  equation  solutions  with  exact 

normal  mode  solutions.  For  a  particular  mode  of  order  m  with  phase 

velocity  c  =  w/k  ,  and  initial  range  rn,  Fitzgerald  suggested  the 
mm  w 

following  limitation  on  the  maximum  distance,  r  -  r^,  over  which  the 
parabolic  equation  is  useful  at  frequency  f: 

r  -  rn  «  c  2cn  /  [  f ( c  -  cn)2]  .  [8] 

The  above  equation  applies  to  sound  waves  in  the  SOFAR  channel  and  to 
RSR  (refracted  surface  reflected)  propagation.  Typical  limits  imposed 
by  [8]  are  111  km  at  100  Hz  to  15,000  km  at  10  Hz.  For  the  case  of 
bottom  limited  propagation,  the  required  small  finite  difference  step 
size  and  consequently,  the  long  time  spent  computing,  poses  a  greater 
practical  restriction  on  the  valid  range  than  equation  [8].  In  general, 
the  parabolic  equation  will  be  accurate  to  longer  ranges  as  the 
frequency  f  decreases  or  c^  approaches  Cq. 
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An  implicit  finite  difference  formulation  of  the  parabolic  equation 


Following  McDaniel  [1975b],  an  implicit  finite  difference  algorithm 
is  constructed  for  equation  [7].  To  simplify  the  notation,  equation  [7] 
can  be  written  as 


u  =  u  +  bu  ,  [  9  ] 

zz  —  r 
a 

where  a  =  i/2k,  b  =  2k(kQ  -  k),  and  k  =  k(z,r)  can  be  complex. 

Now,  equation  [9]  is  solved  numerically.  In  the  discussion  that 
follows,  the  subscript  m  indexes  steps  in  the  z  direction  and  the 
superscript  n  indexes  steps  in  the  r  direction.  A  grid  point  has 
the  coordinates  (m,n)  and  the  value  of  the  wavefield  at  that  point  is 

u11  .  As  n  is  incremented,  the  wavefield  is  propagated  in  the  r 
m 

direction.  At  each  range  step  n,  the  algorithm  moves  through  values  of 
m  in  the  z  direction.  The  first  derivative  in  r  is  approximated  by 


where  d  is  the 
n 

approximate  local 


u 

r 


n+1  n 

u  -  u 
m  m 


0(d  ) 
n 


n 

range  grid  spacing  from  u11  to 
truncation  error  at  point  (m,n). 


n+1 

u 


[10] 

0(d  )  is  the 
n 


Since  a  difference 
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expression  such  as  equation  [10]  is  derived  from  a  Taylor  series,  the 
truncation  error  results  from  the  neglected  terms  in  the  series.  The 

second  derivative  in  z  at  u11  is  approximated  by  [Gerald,  197  0] 


n  n 

u  .  -  u 
0  m+1  m 

u  =*2  _  - 

ZZ  H  h 

m 

where  11=  (h  +  h  .),  h  =  H/2,  and  h 
m  m— 1  m 

u  to  u  .  •  When  every  h  is  the  same 
m  m+1 

[11]  takes  on  the  common  form 


n  n 
u  -  u  , 

2  +  0(h) 

E  h  , 
m-1 

is  the  depth  grid 
(equispaced  grid  in 


[11] 

spacing  from 
z) ,  equation 


u 


zz 


u*  -  2U11  +  Un  .  ? 

m+1  m  m~1  +  o(h2)  . 


[12] 


Note  that  the  error  has  fallen  to  0(h  ). 

For  the  moment,  we  shall  let  h  be  constant.  The  second  difference 

2 

operator  6^  is  defined  as  follows: 


,2  n 
o  u 
z  m 


*  u 


-  2u 


J.  -*  T  U  1  . 

m+1  m  m-1 

The  actual  second  derivative  and  the  second  difference  operator  have  the 


[13] 


following  relation: 


u 

zz 


n 

o  u  0 

z  m  +  0(h2) 


[14] 
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From  a  Taylor  expansion  [Mitchell  and  Griffiths,  1980], 

62un 

z  m  «  (u  )”  +  1  h2(u  )°  +  0(h4)  .  [15] 

— —  zz  m  YJ  zzzz  m 

h 

Averaging  derivatives  between  the  two  range  steps  n,  n+1,  gives  a  form 
of  the  frequently  used  Crank-Nico 1  son  equation  for  solving  finite 
difference  problems : 


z  /  n+1  nN 

_ Cu  +  u  ) 

r,  m  m 

2h 


u  ri/2  - 

zz  m 


1  h2(u  )”+1/2 

|  ^  ZZZZ  IQ 


+  0(h4+  d  ) 


[16] 


where  the  field  value  u  averaged  between  the  two  range  steps  n,  n+1,  is 

n+1  ,  n 
i/o  u  +  u 
n+1/2  m  m 

u  =  _  . 

m  2 

The  left  side  of  equation  [16]  could  more  generally  be  written  as 


£ 

Z  (  eun+1  +  (1  -0  )un)  .  [17] 

— r  m  m 

h 

In  the  case  of  equation  [16],  0=  1/2.  0=1  yields  the  fully  implicit 

formula  while  using  0=0  results  in  the  fully  explicit  formula.  For 
1/2  <  0  <  1,  this  type  of  differencing  scheme  is  unconditionally  stable, 
even  for  the  case  of  complex  coefficients  [Lee  and  Fapadakis,  1979]. 
Hood  [1978],  and  Claerbout  [1970a]  state  that  using  a  value  of  6  greater 
than  1/2  may  improve  the  accuracy  of  the  propagated  wavefield.  However, 
1/2  is  a  practical  choice  since  the  increased  accuracy  obtained  by 
selecting  another  value  is  usually  marginal  when  considering  the 
computer  time  required  to  determine  the  optimum  value. 
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Substituting  the  right  side  of  equation  [9]  as  a  replacement  for 

u  in  the  right  side  of  equation  [16]  and  also  using  equation  [10] 
zz 

produces 
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Now  we  let  h  =  H/ 2  and  use  equation  [11]  to  define  the  6^  operator  for 
a  grid  with  variable  spacing; 
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The  final  difference  equation  is  obtained  by  inserting  equatior 
into  equation  [18]  and  rearranging  to  get  equation  [20]; 
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^  2 

The  local  truncation  error  of  equation  [20]  is  0(h  +  d  )  when  the 

finite  difference  grid  is  equispaced.  Equation  [20]  may  be  rewritten  as 


.n+1  n+1  ,  _n+l  n+1  .  „n+l  n+1 

A  .  u  .  +  B  u  +  C  .  u  ,  = 

m-1  m-1  m  m  m+1  m+l 

Dn  .  un  .  +  En  un  +  F“  un  .  . 

m—  1  m- 1  mm  m+ 1  m+ 1 


[21] 


m  runs  from  1  to  M  where  M  is  the  total  number  of  grid  points  in  the  z 
direction.  The  set  of  equations  represented  by  equation  [21]  may  be 
written  as  the  tridiagonal  matrix  equation  [22]; 


D1  U1  +  E2  u2  +  F3  u3  +  T0P 


ri  _,n  n  _n  n 

D0  u0  +  E0  u0  +  F.  u. 

2  2  3  3  4  4 


-ji  n  n  n  n  n 
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_n  n  ^  _n  n  ^  „n  n 
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_n  n  n  „n  n 

D  qU  «  +  E  0u  0+  F  -a  u  , 
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D  0u  0+  E  iu  i+  F  u  + 
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The  wave  field  is  propagated  along  the  r  direction.  For  each  range 
index  n,  the  matrix  or  equation  [22]  must  be  inverted  to  solve  for  the 
field  u  at  index  n+1  •  The  term  "implicit  finite  difference"  comes  from 

the  fact  that  the  field  value  un+^  is  determined  from  the  field  at  n  by 
inverting  a  matrix  rather  than  from  direct  multiplication.  Special 
numerical  methods  are  available  to  invert  tridiagonal  matrices  such  as 
the  left  side  of  equation  1 22 J  [Hornbeck,  1975]. 

Since  m  =  1  corresponds  to  the  surface  of  the  ocean,  the  free 
surface  (pressure  equal  to  zero)  boundary  condition  is  implemented  by 

setting  Uj  =  0  and  the  term  TOP  =  0  for  all  n.  A  simplified  version  of 
the  Neumann  boundary  condition  given  by  Lee  and  Papadakis  [1979]  is  used 
at  the  very  bottom  of  the  grid.  Note  that,  in  a  homogenous  medium  with 
a  uniformly  spaced  grid  and  zero  pressure  at  both  surface  and  bottom, 
the  left  hand  matrix  of  equation  [22]  is  symmetric. 

The  accuracy  of  the  implicit  finite  difference  algorithm  developed 
here  was  verified  by  computing  the  steady  state  wave  field  for  published 
models  [Volk,  1975;  Jensen  and  Kuperman,  1980J. 
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Source  modeling 

The  acoustic  data  presented  in  the  next  chapter  consist  of  many 
shots  recorded  by  one  receiver.  To  model  such  data,  the  acoustic  source 
is  placed  at  the  location  of  the  original  receiver.  The  principle  of 
reciprocity  states  that  propagation  from  the  position  of  the  receiver  to 
each  shot  is  the  same  as  that  for  shot  to  receiver.  The  source  used  in 
the  modeling  is  the  Gaussian  pulse  specified  by  Tappert  [1977]; 

Source(z,r=0)  =  P  exp(  -(z  -  zg)^  /  (2/1Cq))  •  [23] 

P  =  pQ  /  (  /2/k Q) 

zg  =  source  depth 
Pq  *  initial  pressure 

The  exact  value  of  the  source  amplitude  P  does  not  matter  since  absolute 
signal  levels  are  not  modeled. 


Including  attenuation  in  the  parabolic  equation 

Attenuation  is  introduced  by  using  a  complex  wavenumber,  (k  +  iot )  , 
where  k  *w/c.  The  standing  wave  exp(ikz)  along  the  z  axis  becomes, 
with  attenuation; 


i(k  +  ia)z  -az  ikz 

e  =  e  e 


[24] 
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The  units  of  a  are  nepers  (Np)  per  unit  length  [Sheriff,  1973],  If  A 
is  the  attenuation  coefficient  in  decibels  per  unit  length,  then, 


A  A 

a  =  _  ■  _ 

20  log^e  8.686 

[Clay  and  Medwin,  1977].  In  this  thesis,  it  is  assumed  that 
linearly  with  frequency  [Hamilton,  1972], 

Attenuation  expressed  in  terms  of  the  intrinsic  Q  is 
Richards,  1980] 


[25] 

A  varies 

[Aki  and 


exp(  ^  )  =  exp(  -ctz) 

2  c  Q 

1  -  If.  [26] 

Q  f 

where  c  is  the  phase  velocity  and  f  is  frequency*  For  a  25  Hz  wave 
traveling  at  1500  m/ s ,  an  a tt enua tion  coef f icient  of  1.5  dB/km  is 
equivalent  to  a  Q  of  about  100. 


Modeling  the  ocean-sediment  interface 

The  modeling  of  discontinuities  such  as  the  water-sediment 
interface  requires  special  attention.  Both  the  density  and  velocity 
contrasts  across  a  boundary  must  be  considered.  The  use  of  the 
acoustic  wave  equation  [1]  requires  that  the  ocean  bottom  be  modeled  as 
fluid  rather  than  solid  material.  The  assumption  of  a  fluid  bottom  is 
thought  to  be  reasonable  for  many  ocean  acoustics  problems,  except  when 
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the  sediments  are  thin  (Vidmar,  1980)*  In  thin  sediments,  the  P  to  S- 
wave  conversions  at  the  sediment-basement  interface  cause  propagation 
loss  effects  that  are  ignored  by  a  fluid  sediment  model. 

McDaniel  and  Lee  [1982]  developed  an  implicit  finite  difference 
scheme  that  specifically  satisfies  acoustic  interface  boundary 
conditions,  but  the  method  is  not  readily  adaptable  to  the  case  of  a 
sloping  bottom.  Interfaces  can  be  approximated  by  rapid,  but 
continuous,  variations  in  velocity  and  density,  and  that  is  the  method 
used  here.  Replacing  the  wavenumber  k  (k  =  k(r,z))  of  equation  [7]  with 
K  given  below  will  include  the  effect  of  density  gradients 
[Brekhovskikh,  1980,  p.  162;  Tappert,  1977]. 

2  ,  2  1  —.2  3  / 1  %  2  r  07 1 

K  =  k  +  _ V  p  -  _  (_Vp)  [27] 

2  P  4  P 

The  terms  of  equation  [27]  involving  density  may  be  neglected  if  the 
changes  in  density  are  small  over  a  wavelength.  One  difficulty  in  using 
equation  [27]  is  that  K  is  very  dependent  upon  the  curvature  of  the 
rather  arbitrary  curve  that  one  chooses  to  represent  the  density 
contrast  at  an  interface.  Rather  than  explicitly  varying  k  through  the 
use  of  equation  [27],  k  (or  c)  in  the  vicinity  of  an  interface  will  be 
specified  so  that  the  theoretical  reflectivity  of  the  interface  region 
is  close  to  the  Rayleigh  reflectivity  of  the  true  acoustic  interface. 

The  theory  of  Brekhovskikh  [1980]  was  implemented  to  model  an 
interface  by  simultaneously  considering  the  density  and  velocity 
contrast  across  the  interface.  The  interface  is  approximated  by  an 
Epstein  layer  (named  after  the  Soviet  physicist  P.  Epstein)  with 
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2 

parameters  M,  N,  and  m.  The  square  of  the  index  of  refraction,  n  (2), 
for  a  horizontal  Epstein  layer  where  k  -  k(z)  and  c  -  c(z)  is  given  by 


2 
n 

N  exp(  m(z  -  zn))  4  M  exp(  m(z  -  zn)) 

1  -  _ °  -  _ _  .  [28] 

2 

(1  +  exp(  m(z  -  Zq)  )  (1  +  exp(  m(z  -  z^)) 

Zq  is  the  Epstein  layer  depth  (z  measured  positive  downward),  is  the 
constant  velocity  in  the  medium  well  above  the  layer,  and  k^  =  2  tt  /  X  ^  . 
The  constant  M  depends  upon  the  density  contrast  at  the  interface,  while 
N  is  determined  from  the  velocity  contrast.  A  plane  wave  propagates  in 
the  medium  above  the  layer  and  is  both  reflected  and  transmitted.  In 
the  present  work,  acoustic  interfaces  are  modeled  using  equation  [28]. 
Setting  M  =  0  gives 


n  N  exp(  m(z  -  zn)) 

n2  =  1  -  _ °  ,  [29] 

(1  +  exp(  m(z  -  z^)) 


which  represents  the  square  of  the  index  of  refraction  due  solely  to  a 
velocity  contrast  across  the  interface  (Figure  1  ).  For  z  <<  z^ 

2  2 

equation  [28]  becomes  n  =  (c^  /  c)  =  1 ,  i.e.,  c(z)  =  c^  above  the 

2  2 

Epstein  layer.  For  z  »  Zq  n  =  (c^  /  c^)  *  1  -  N  ,  where  c^  is  the 
velocity  below  the  Epstein  layer.  Thus,  N  is  determined  from  N  *  1  - 


(cx  /  c2r 


Fig,  1,  Left:  n  ,  the  index  of  refraction  squared,  for  the  symmetrical 
component  (equation  [30])  of  the  Epstein  layer  which  describes  the 
density  contrast  across  an  interface  (M  =  0,4), 

2 

Right:  n  for  the  transitional  component  (equation  [29])  of  the  Epstein 

layer  which  accounts  for  the  velocity  contrast  across  an  interface 
(c^  =  1490  m/s,  =  1641  m/s).  Parameters  common  to  both  figures  are 

l  =  10  m,  and  zn  =  50  m. 
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Now,  setting  N  «  0  in  equation  [28]  results  in  the  equation  for  a 
density  contrast  only 


2 

n 


1 


4  M  exp(  m(  z  -  Zq) ) 

-- _  * 


[30] 


2 

( 1  +  exp(  m( z  -  Zq) ) 

which  is  plotted  in  Figure  1.  Note  that  equation  [30]  is  near  unity 
except  when  z  is  near  Zq  ,  where  it  approaches  1  -  M.  Morris  et  al. 
[1978]  suggested  iteratively  changing  the  model  velocity-depth  function 
near  the  continuous  interface  until  the  model  reflectivity  matches  the 
Rayleigh  reflectivity  of  the  true  interface*  Here,  modeling  interfaces 
proceeds  along  similar  lines.  In  this  case,  M  is  varied  by  trial  and 
error  until  the  Epstein  layer  reflectivity  matches  the  Rayleigh 
reflectivity  for  the  interface.  Thus,  M  is  not  directly  specified  from 
the  known  density  values  above  and  below  an  interface. 

The  term  m  is  related  to  the  Epstein  layer  thickness  and  may  be 
interpreted  using  equation  [30].  The  distance  between  the  half  maximum 
points  of  equation  [28]  is  ^  =  3.5254/m  •  The  parameter  S  =  2k^/m  is 
called  the  relative  thickness  of  the  layer,  k^  «  2  /X  where  X ^  is  the 
incident  wavelength.  A  should  be  smaller  than  X ^ ,  yet  larger  than  the 
computational  grid  spacing. 

Based  upon  a  solution  of  the  reduced  wave  equation  [2], 

Brekhovskikh  [1980,  p.  171]  computed  the  complex  reflection  coefficient 

or  R  (P-wave  incident,  P-wave  reflected),  for  the  Epstein  layer  of 
PP 

equation  [28]; 
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R  =  r(y- l)  r(i -s  )  r  (i  +  a -y  )  ^  [31] 
r(i  -y  )  r(  y  -  8  )  r  (a) 

a,  8,  and  y  are  complex  functions  of  the  angle  of  incidence,  wavelength 
of  incident  wave,  m,  N,  and  M  and  are  given  in  Appendix  A,  The  T  (  ) 
are  complex  gamma  functions  which  are  numerically  evaluated  by  use  of  a 
seven— term  series  expansion  due  to  Lanczos  [Luke,  1969]#  A  Fortran 
program  to  compute  the  complex  gamma  function  is  provided  in  Appendix  A. 

The  Rayleigh  reflection  coefficient  for  an  acoustic  interface  where 
the  subscript  1  indicates  a  property  of  the  upper  medium,  2  refers  to 
the  lower  medium,  and  ®  is  the  angle  of  incidence,  is 


R 

PP 


X  -  Y 
X  +  Y 


p2 

where  X  =  _  ,  and  Y 

pi 


(  (c1  /  c2)2  -  sin2©)1^2 


•  2,1/2 

(1  -  sin  6) 


[32] 


The  interface  reflectivity  computed  from  equation  [32]  provides  the 
model  to  which  the  Epstein  layer  reflectivity  is  matched. 
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III.  MODELING  PROPAGATION  ON  A  CONTINENTAL  SLOPE. 


Field  experiment 

The  acoustic  signals  modeled  in  this  portion  of  the  study  were 
recorded  near  the  edge  of  the  continental  shelf  off  Nova  Scotia  in  June 
1975.  All  signals  were  received  by  ocean  bottom  seismometers,  each 
equipped  with  a  triaxial  set  of  geophones  and  a  hydrophone.  The 
instrument  package,  known  as  a  TOBS  (Telemetering  Ocean  Bottom 
Seismometer),  has  been  described  by  Sutton  et  al.  [1977].  In  the 
succeeding  sections,  propagation  loss  modeling  results  for  SUS  (sound 
underwater  signaling)  explosive  charges  and  airgun  shots  will  be 
discussed.  The  hydrophone  data  was  considered  here  since  the  parabolic 
equation  calculates  the  acoustic  pressure  field  for  an  environmental 
model.  The  hydrophone  frequency  response  was  flat  from  about  10  Hz  to 
40  Hz.  Because  acoustic  pressure  and  particle  velocity  are  linearly 
related  [Clay  and  Medwin,  1976,  p.  59],  rates  of  propagation  loss  with 
range  for  acoustic  signals  measured  by  geophones  and  hydrophones  should 
generally  be  identical. 

The  SUS  charge  line  was  situated  downslope  from  the  ocean  bottom 
seismometer  (TOBS  3)  on  an  average  slope  of  3.5°  (Figure  2).  At  TOBS  3 
the  water  depth  was  1301  m,  while  at  the  furthest  shot,  the  water  depth 
was  over  2800  m.  The  shots  were  1.1  oz  (0.031kg)  SUS  charge  boosters 
detonated  at  a  depth  of  18  m.  Figure  2  shows  the  location  of  the  SUS 


Fig.  2.  Location  of  TOBS  3  and  SUS  charge  line  on  the  Scotian 
continental  slope.  Velocimeter  stations  4  and  5  are  the  sources  of 
water  velocity  data  for  the  SUS  charge  line.  Also  shown  is  the  airgun 
line  1-112  and  receiver  TOBS  1  on  the  shelf. 
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charge  line  on  the  continental  slope,  as  well  as  the  position  of  two 
sound  velocity  profiles  (Figure  3), 

The  SUS  charge  arrivals  are  plotted  in  Figure  4.  Also  shown  in  the 
figure  is  the  6.4  second  data  window  from  which  the  propagation  loss 
rates  were  calculated.  The  signal  strength  was  greatest  near  25  Hz, 
which  corresponded  to  the  bubble  pulse  frequency.  Individual  bottom 
bounce  multiples  are  easily  identified  in  Figure  4.  The  highest 
discernable  group  velocity  was  close  to  1500  m/s,  indicating  that  most 
of  the  energy  propagated  in  the  water  column. 


Physical  properties  of  the  slope 

In  general,  acoustic  modeling  requires  detailed  environmental 
input  because  it  is  well  known  that  the  properties  of  the  uppermost 
sediments  greatly  influence  acoustic  wave  propagation  in  shallow  water. 
Unfortunately,  little  sediment  information  for  this  SUS  charge  line  was 
available.  A  1.8  lb  (0.82  kg)  SUS  charge  refraction  line  over  the  same 
location  revealed  an  upper  sediment  velocity  of  about  1.7  km/ s  [Brocher, 
1982,  in  prep.].  Ray  tracing  [Appendix  B]  identified  a  bottom  incident 
critical  angle  of  about  65°.  This  result  was  supported  by  the 
observation  of  a  head  wave  at  close  range  resulting  from  incidence  at 
65.2°  [Appendix  C].  A  bottom  water  sound  speed  of  1490  m/ s  near  TOBS  3 
and  critical  angle  of  65.2°  yielded  an  upper  sediment  velocity  of 
1641  m/s. 
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Fig*  4.  SUS  charge  arrivals  on  the  Scotian  continental  slope  from  shots 
downslope  of  the  receiver.  The  dashed  lines  show  the  6.4  second  data 
window  used  for  propagation  loss  measurement.  In  addition  to  the  direct 
wave,  three  bottom  multiple  arrivals  are  observable. 
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Piston  cores  up  to  13  m  long  from  the  Scotian  continental  slope 
contained  muds  and  clays  with  a  small  percentage  of  sand  [Piper,  1974]  . 
The  1641  m/s  mentioned  above  was  closer  to  the  compressional  velocity  in 
silt  rather  than  to  the  lower  velocity  in  clay  [Hamilton,  1980].  The 
total  thickness  of  the  unconsolidated  sediments  was  not  determined.  On 
the  nearby  shelf,  the  unconsolidated  sediments  were  on  the  order  of  tens 
of  meters  thick,  occasionally  thickening  to  over  200  meters  [Parrott  et 
al.,  1980] . 


Constructing  a  velocity-depth  model 

A  velocity-depth  function  was  obtained  from  the  velocimeter 
information  shown  in  Figure  3  and  upper  sediment  velocities  determined 
from  the  seismic  refraction  data.  The  water  velocity-depth  function 
used  in  the  geoacoustic  model  was  a  range  weighted  average  of  data  from 
sound  velocity  stations  4  and  5  between  0  and  4  km  range,  and  was 
identical  to  the  station  5  sound  velocity  profile  beyond  4  km. 

An  ocean  bottom  sediment  velocity  of  1641  m/s  was  selected,  based 

on  the  critical  angle  estimate.  A  sonic  velocity  gradient  of  1  s  ^  was 
chosen  for  the  sediments.  This  gradient  lies  within  the  range  of  values 
reported  for  silt  clays  and  turbidites  [Hamilton,  1979].  The  sediment 
thickness  was  arbitrarily  set  at  60  m.  As  a  result  of  using  this 
sediment  velocity  and  thickness,  and  a  bottom  water  velocity  of  1492  m/s 
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(measured  about  7  km  downslope  of  TOBS  3),  the  minimum  bottom  incidence 
angle  for  long  range  propagation  was  61.3°.  That  is,  rays  with  a  bottom 
incidence  of  61.3°  or  more  remained  in  the  water  column  either  by 
reaching  a  turning  point  within  the  upper  60  m  of  sediment  or  by 
postcritical  reflection  at  the  sediment-water  interface.  In  the  model, 
bottom  incident  rays  at  angles  less  than  61,3°  were  unable  to  turn  back 
upwards  within  the  upper  60  m  of  sediment  and  were  lost  from  the  sound 
channel • 

A  special  provision  was  made  for  rays  refracted  from  the  sound 
channel  into  the  bottom  of  the  sediment.  Although  it  might  be 
preferable  to  follow  Clayton  and  Engquist  [  1977],  in  which  absorbing 
boundary  conditions  were  used  to  allow  waves  to  pass  out  of  the 
computational  grid  without  any  reflection  or  backscatter,  it  is 
nontrivial  to  apply  the  Clayton  and  Engquist  boundary  condition  on  a 
slope*  In  this  case  the  waves  passing  through  the  bottom  of  the 
sediments  were  simply  damped  out  as  was  done  by  Tappert  [1977]  *  More 
specifically,  for  every  other  depth  grid  point  below  the  sediment  layer, 
the  attenuation  was  doubled  while  the  velocity  was  remained  constant. 

This  doubling  continued  until  the  final  attenuation  was  2^  times  the 
attenuation  in  the  sediments.  Since  the  wavefield  was  damped  out  by  the 
time  it  reached  the  bottom  of  the  computational  grid,  the  boundary 
condition  at  the  very  bottom  did  not  matter,  as  long  as  it  was 
physically  realizable.  The  Neuman  boundary  condition  on  a  slope  was 
specified. 
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A  sediment  density  of  1.6  g/cm  and  an  attenuation  in  the  sediments 
of  1.5  dB/km  at  25  Hz  were  selected  for  use  in  the  model  from  values 
computed  by  Beebe  and  McDaniel  [1980]  for  the  Scotian  shelf.  Beebe  and 
McDaniel  used  sediment  grain  size  data  and  the  Biot-Stoll  sediment  model 
[Stoll,  1974]  to  compute  the  attenuation.  These  parameters  are 

approximate,  since  it  was  likely  that  the  Scotian  shelf  sediments  had  a 
higher  sand  content  than  those  on  the  Scotian  slope.  Figure  5 
summarizes  the  environmental  parameters  used  in  the  model. 


An  Epstein  layer  for  the  water-sediment  interface 

The  ability  of  ray  tracing  in  a  medium  with  a  hard  bottom  to 
qualitatively  predict  variations  in  signal  level  [Appendix  B]  suggested 
that  the  water-sediment  interface  was  the  principal  bottom  reflector  and 
consequently  a  deeper  acoustic  basement,  if  any,  did  not  significantly 
affect  the  guided  wave  propagation  in  the  water.  Therefore,  the  only 
bottom  reflector  included  in  this  parabolic  equation  model  was  the  ocean 
bottom.  As  described  earlier,  the  ocean  bottom  was  simulated  by  an 
Epstein  layer. 

Figure  5  shows  the  geoacoustic  model  and  the  Epstein  layer 
approximation  to  the  ocean  bottom  interface.  An  Epstein  layer  half¬ 
thickness  A  ®  10  m  was  chosen.  This  choice  was  subject  to  the 
constraints  that  it  be  greater  than  the  computational  grid  size  of  2  m 
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Fig,  5.  Left:  Geoacoustic  parameters  for  modeling  propagation  on  the 
Scotian  slope. 

Right:  The  equivalent  Epstein  layer  model  of  the  ocean  bottom 
interface* 
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in  the  z  direction,  and  less  than  the  average  wavelength  of  60.8  m 
computed  from  a  25  Hz  wave  traveling  with  an  average  group  velocity  of 
1520  m/s.  Once  i.  (or  small  m)  was  selected,  this  left  the  parameter 
M,  associated  with  the  density  contrast,  as  the  only  free  variable  for 
matching  the  Epstein  layer  reflectivity  [32]  to  the  Rayleigh 
reflectivity  [34].  The  Rayleigh  and  Epstein  reflectivities  calculated 
from  the  parameters  shown  in  Figure  5  are  compared  in  Figure  6.  The 
choice  M  =  0.4  was  based  upon  trial  and  error  matching  of  the  magnitudes 
of  the  respective  reflection  coefficients.  Note  that  while  the 
magnitude  of  the  reflection  coefficient  was  reasonably  well  matched, 
there  was  a  considerable  mismatch  in  phase  (Figure  6). 
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Fig.  6.  Rayleigh  reflection  coefficient  and  phase  (a)  for  the  ocean 
bottom  of  Figure  5.  Solid  line  shows  the  25  Hz  reflectivity  for  the 
Epstein  layer  model  of  the  bottom  interface* 
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Choosing  a  computational  grid 

Brock  [1978]  described  a  numerical  algorithm  for  the  parabolic 
equation  which  automatically  picked  the  range  increment  and  z  increment 
as  the  computation  proceeded.  However,  the  algorithm  was  based  upon  a 
Fourier  transform  solution  of  the  wave  equation  and  uniform  z  increment 
at  each  range  step  and  was  not  directly  applicable  to  this  study.  In 
the  present  analysis,  the  computational  grid  was  arbitrarily  chosen 
using  the  following  criteria: 

i)  There  must  be  several  grid  points  per  wavelength  of  the 
propagating  waves. 

ii)  The  grid  increment  must  be  smaller  than  the  "wavelength”  of  the 
boundaries  as  well  as  of  the  source. 

Both  of  the  above  requirements  are  restatements  of  the  fact  that  the 
grid  step  size  must  be  small  enough  to  prevent  spatial  aliasing. 

Table  1  lists  the  specifications  for  a  fixed  and  variable  step  size 
grid  on  which  the  model  described  in  Figure  5  was  computed.  The  results 
for  both  grids  are  compared  in  the  next  section.  The  relative  grid 
spacings  are  illustrated  in  Figure  7. 

It  is  necessary  to  specify  a  small  depth  increment  near  the  ocean 
bottom.  In  a  rectangular  grid,  this  small  increment  in  depth  persists 
to  all  ranges.  With  a  sloping  bottom,  the  result  is  a  fine  grid  where 
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Table  1. 


—  Finite  difference  grid  spacings  for  modeling  propagation  on 
the  continental  slope  near  Nova  Scotia. 


Grid  Spacing  in  the  Depth  Direction 


Depth, 

Grid  Spacing 

Index 

meters 

Az.  meters 

i 

0 

2 

50 

98 

51 

100 

4 

320 

1176 

321 

1180 

2 

1200 

2938 

Grid  Spacing 

in  the  Range 

Direction 

Range, 

Grid  spacing 

Index 

meters 

Ar.  meters 

1 

0 

5 

1200 

5995 

1201 

6000 

10 

2400 

17990 

2401 

18000 

20 

2919 

28360 

Az  =  2  meters  and  Ar  =  5  meters  for  the 
constant  spacing  grid. 


Fig,  7.  Above:  Constant  grid  spacing  model. 

Below:  Variable  grid  spacing.  For  both  grids  only  about  1  in  40  grid 
intersections  in  both  range  and  depth  direction  are  plotted.  Note  that 
the  grid  length  in  the  z  direction  varies  with  the  bathymetry.  The 
complete  grid  specifications  are  listed  in  Table  1. 
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it  is  not  needed  and  increased  computational  costs.  There  are 
alternatives  to  the  rectangular  variable-spacing  grid  used  in  this 
thesis.  If  the  slope  is  fairly  uniform,  as  in  the  present  case, 
effecting  a  non-rectangular  grid  which  conforms  to  the  average  bottom 
slope  saves  computer  time  without  loss  of  accuracy.  Another 
possibility  is  to  change  the  finite  difference  grid  at  a  specified 
range.  At  that  range,  propagation  is  halted  and  the  wave  field  is 
interpolated  onto  a  new  z  grid.  Computing  the  wave  field  would  then 
proceed  until  it  is  necessary  to  change  the  grid  again. 
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The  parabolic  equation  model  output 

Figure  8  shows  the  intensity  of  the  acoustic  wave  field  for  25  Hz 
generated  by  the  parabolic  equation  model  based  on  Figure  5  and  the 
fixed  grid  spacing  detailed  in  Table  1 •  A  Cq  value  of  1520  m/ s  was 
used  in  the  computation.  In  both  the  range  and  depth  directions,  the 
values  of  the  wavefield  at  about  1  in  16  computational  grid  points  were 
plotted.  At  short  ranges,  the  acoustic  intensity  was  highly  oscillatory 
due  to  the  propagation  of  many  modes  associated  with  different  wave 
numbers  (or  differing  phase  velocities).  With  increasing  distance, 
precritical  reflections  were  quickly  damped  out,  leaving  the 
postcritical  phases  indicated  by  the  ray  traces  of  Figures  17-19. 

Figure  9  (bottom)  is  a  plot  of  the  wavefield  of  Figure  8  at  18  m 
depth,  the  same  depth  as  the  SUS  charge  detonation.  A  least  square  fit 
straight  line  from  6  to  28  km  yielded  a  propagation  loss  rate  (negative 
of  the  slope)  of  0.61  dB/km.  This  value  was  automatically  corrected  for 
cylindrical  spreading  because  of  the  assumption  (equation  [3])  made 
during  the  derivation  of  the  parabolic  approximation.  The  propagation 
loss  for  the  same  model  computed  using  the  variable  grid  spacing  of 
Table  I  was  0.60  dB/km  (  Figure  10).  Figure  10  compares  the  variable  and 
fixed  grid  computations  at  18  m  depth.  A  close  examination  of  Figure  10 
revealed  only  minor  differences,  which  tended  to  increase  with  range, 
between  the  outputs  of  the  two  grids.  The  variable  grid  spacing  was 
more  economical  since  computing  the  same  wavefield  as  in  Figure  8  took 
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Fig.  8.  Intensity  of  25  Hz  parabolic  equation  wave  field.  Input  is  the  model 
of  Figure  5  computed  on  the  constant  spacing  grid.  White  is  the  weakest 
amplitude  and  black  is  the  strongest.  From  black  to  white  inclusive  there  are 
seven  shades  each  separated  by  2.6  dB.  The  Gaussian  pulse  source  is  seen  at  far 
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Fig.  9.  Top:  Observed  25  Hz  SUS  charge  signal  level  versus  range 
computed  using  the  6.4  second  data  window  shown  in  Figure  4.  The 
decibel  values  have  been  corrected  for  cylindrical  spreading.  The 
-0.19  dB/km  propagation  loss  or  negative  slope  is  computed  for  ranges 

from  4  to  23  km. 

Bottom:  Computed  parabolic  equation  wavefield  at  25  Hz  using  the 
constant  grid  spacing  of  Table  1*  This  plot  is  a  slice  taken  at  18  m 
depth  from  Figure  8.  The  +0.61  dB/km  propagation  loss  is  measured  for 
ranges  from  6  to  28  km. 
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25  HZ 

PROP.  LOSS 


CONSTANT  SPACING  A  +0.60  dB/km 
VARIABLE  SPACING  X  +0.61 


RANGE  KM 


Fig.  10.  Comparison  of  the  outputs  of  the  parabolic  equation  model  at 
25  Hz  for  the  variable  and  constant  grid  spacing  listed  in  Table  1  and 
illustrated  in  Figure  7.  The  measured  +0.60  dB/km  propagation  loss  from 
6  to  28  km  computed  using  variable  grid  spacing  is  almost  the  same  as 
the  +0.61  dB/km  loss  for  the  constant  grid  spacing  case* 
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only  33  minutes  of  computing  on  a  Harris  800  computer  and  required 
smaller  arrays  than  the  fixed  grid  calculation  which  lasted  96  minutes. 


Comparison  of  the  observed  and  predicted  propagation  loss 

The  predicted  0.61  dB/km  propagation  loss  at  25  Hz  does  not  match 
the  observed  negative  propagation  loss  of  -0.19  dB/km  (Figure  9).  It 
was  not  known  why  the  signal  level  increased  with  range  after 
normalizing  for  cylindrical  spreading  loss.  A  possible  problem  was  that 
the  instrumentation  compressed  the  hydrophone  signals,  so  that  the 
apparent  propagation  loss  was  less  than  expected  for  cylindrical 
spreading. 

Comparing  the  observed  signal  levels  in  Figure  9  with  the  parabolic 
equation  predictions  in  the  same  figure  showed  possible  common  peaks  at 
6.5,  9.5,  16,  and  23  km.  In  general,  the  shot  spacing  was  too  coarse  to 
enable  resolution  of  sudden  changes  in  signal  level.  Significant 
arrivals  predicted  at  11.5  to  12.5  km  and  16 .5  to  18.5  km  were  absent 
from  the  data.  Errors  in  the  geoacoustic  model  and  in  modeling  the 
source,  as  well  as  the  signal  compression,  may  explain  the  discrepancy 


between  data  and  model. 
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IV.  MODELING  DATA  FROM  A  CONTINENTAL  SHELF 


Bathymetry 

Airgun  profiling  by  an  oil  company  vessel  provided  signals  along  a 
line  (1-112)  oblique  to  TOBS  1  (Figure  2).  This  data  set  provided  the 
opportunity  to  analyze  sound  propagation  in  shallow  water  over  a  slope 
of  less  than  0.1  degrees.  In  addition,  the  data  were  of  good  quality 
and  the  shot  coverage  was  dense  enough  to  resolve  sharp  signal 
fluctuations  over  distances  of  a  few  hundred  meters.  The  bathymetry 
along  the  airgun  line  is  plotted  in  Figure  11.  Because  of  the  rapidly 
increasing  water  depth  past  29*5  km  only  data  at  ranges  less  than  29*5 
km  were  used*  The  water  depth  varied  from  about  90  m  at  31.6  km  to  83  m 
at  a  range  ot  10.3  km.  As  Figure  2  indicates,  the  bathymetry  between 
the  shot  locations  and  TOBS  1  was  three  dimensional. 

To  model  the  bathymetry,  the  ocean  bottom  for  ranges  between  0  and 
25.4  km  was  approximated  by  a  linear  increase  in  depth  from  67  to  70  m. 
For  ranges  from  20  km  to  the  nearest  shot  at  10.3  km,  this  model 
bathymetry  was  an  inaccurate  representation  of  the  actual  water  depth 
beneath  each  shot.  The  largest  error  was  a  difference  of  15  m  at  10.3 
km  range.  Nonetheless,  the  bathymetry  model  provided  a  good 
representation  of  the  water  depth  for  travel  paths  between  the  receiver 
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Fig.  11.  Bathymetry  beneath  airgun  shot  points. 
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Constructing  a  velocity-depth  model 

The  unconsolidated  sediments  in  the  vicinity  of  Banquereau  Bank 
are  glacial  and  littoral  deposits  [Maclean  and  King,  1971;  Parrott  et 
al.,  1980].  The  surficial  unit,  Sable  Island  Sand  and  Gravel,  is 
generally  no  more  than  20  m  thick  [Parrot  et  al.,  1980]  .  Over  much  of 
the  shelf  the  Sable  Island  Sand  and  Gravel  is  underlain  by  the  Emerald 
Silt  and  a  glacial  till  known  as  Scotian  Shelf  Drift.  Both  of  these 
units  vary  in  thickness  from  a  few  meters  to  over  a  hundred  meters. 
Bedrock  is  mostly  Tertiary  sandstones  and  shales  called  the  Banquereau 
Formation  [Jansa  and  Wade,  1974]. 

The  recorded  airgun  signals  are  shown  in  Figure  12.  Most  of  the 
energy  in  the  arrivals  was  concentrated  between  10  and  40  Hz  -  the 
peak  frequency  response  of  the  hydrophone.  A  calculation  based  on  the 
apparent  velocities  of  the  direct  and  trailing  edge  of  the  water  wave 
arrivals  [Houtz,  1980]  [Appendix  D]  yielded  an  upper  sediment  velocity 
of  1.53  km/s.  Measurements  by  McKay  and  McKay  [1982]  using  a  deep-towed 
device  on  the  adjacent  Sable  Island  Bank  found  upper  sediment  velocities 
from  1.57  to  1.66  km/s.  The  sediment  velocity  used  in  the  model  was  1.6 
km/  s. 

Using  a  Biot-Stoll  sediment  model  [Stoll,  1974],  Beebe  and  McDaniel 
[1980]  calculated  a  sediment  attenuation  of  1.5  dB/km  at  25  Hz  (Q=100) 
for  various  locations  on  the  Scotian  shelf.  If  the  exponential 
attenuation  coefficient  varies  linearly  with  frequency  then  the  10  Hz 
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attenuation  is  0.6  dB/km.  The  bottom  density  of  1.8  g/cm  shown  in 
Figure  13  was  taken  from  the  silty  sand  value  given  by  Hamilton  [1980]. 
The  acoustic  basement  on  the  shelf  was  modeled  using  the  parameters  of 
Beebe  and  McDaniel  [1980]  which  were  a  velocity  of  2.00  km/s  and  a 
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density  of  2.2  g/cm  . 

The  acoustic  model  for  propagation  on  the  shelf  contained  two 
interfaces  -  one  at  the  ocean  bottom  and  one  at  the  acoustic  basement. 
The  model  parameters  are  illustrated  in  Figure  13.  Also  shown  in  this 
figure  are  Epstein  layer  models  of  these  interfaces  at  25  and  10  Hz. 
Based  on  an  analysis  of  surface  wave  dispersion  [Brocher,  1982,  in 
prep.],  an  unconsolidated  sediment  thickness  of  21  m  was  used. 

The  25  Hz  reflectivity  functions  of  the  Epstein  layer  equivalents 
to  the  ocean  bottom  and  acoustic  basement  interfaces  of  Figure  13  are 
plotted  in  Figure  14.  As  discussed  earlier,  the  parameters  of  the 
Epstein  layer  were  chosen  by  a  trial  and  error  matching  of  the  Epstein 
reflectivity  curve  to  the  Rayleigh  reflectivity  curve  calculated  from 
the  velocity  and  density  contrast  across  the  interface.  The  10  Hz 
Epstein  layer  reflectivity  plots  for  the  same  interfaces  are  shown  in 
Figure  14  (lower).  At  10  Hz  it  was  necessary  to  use  a  10  m  Epstein 
layer  half  “thickness  compared  to  5  m  at  25  Hz  in  order  to  adequately 
represent  the  acoustic  interfaces.  Thus,  the  extent  of  each  10  Hz 
Epstein  layer  was  20  m  which  is  about  the  same  as  the  model  sediment 
thickness.  Since  the  interfaces  at  the  top  and  bottom  of  the  sediment 
were  modeled  at  10  Hz  by  Epstein  layers  which  were  as  thick  as  the  21  m 


Fig.  13.  Geoacoustic  model  for  the  airgun  shot  line  on  the  Scotian 
shelf.  There  is  an  acoustic  interface  on  the  ocean  bottom  and  at  the 
bottom  of  the  sediment.  Waves  propagating  below  the  subbottom  interface 

are  completely  damped  out. 

Right:  Epstein  layer  approximations  to  the  two  interfaces  at 

25  and  10  Hz. 
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Fig.  14.  Above:  Rayleigh  (A)  and  Epstein  layer  (solid  line) 
reflectivities  at  25  Hz  for  the  ocean  bottom  (upper)  and  sub-bottom 
(lower)  interfaces.  The  left  column  shows  the  magnitude  of  R  ,  the 

right  column  shows  the  phase.  ^ 


Below:  Same  for  10  Hz 
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sediment  layer,  the  10  Hz  modeling  results  must  be  taken  with  caution. 

A  sediment  velocity  gradient  was  not  used  because  of  the  thiness 
of  the  unconsolidated  sediment  layer.  The  water  velocity  was  taken  from 
velocimeter  data  near  TOBS  1.  For  the  computations  a  c^  value  of  1.46 
km/s  was  used.  The  finite  difference  grid  parameters  are  given  in  Table 
2.  The  range  increments  for  computing  the  parabolic  equation  wavefield 
were  the  same  as  for  the  variable  grid  spacing  of  Table  1.  Smaller 
depth  increments  were  used  for  this  model  because  of  the  shallowness  of 
the  water  and  the  need  for  both  bottom  and  sub-bottom  interfaces.  As 
before,  waves  propagating  below  the  sediment  layer  were  damped  out  by 
increasing  the  attenuation  with  depth  while  holding  the  velocity 
constant. 


Comparison  of  data  to  model  at  25  Hz 

The  signal  levels  at  25  Hz  calculated  from  the  6.4  second  data 
window  of  Figure  12  are  plotted  in  Figure  15  (top).  The  observed 
propagation  loss  rate  from  15.4  to  29.5  km  was  -0.17  dB/km.  For  ranges 
less  than  25  km  the  signal  had  a  higher  variance  which  is  probably 
caused  by  the  presence  of  a  number  of  propagating  modes  at  short  range. 
Because  of  the  snallow  water,  most  of  these  modes  rapidly  decayed.  Only 
the  first  mode  could  propagate  beyond  25  km. 
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Table  2 


.  —  Finite  difference  grid  spacings  for  modeling  propagation  on 

the  Scotian  shelf. 


Grid  Spacing 

in  the  Depth  Direction 

Depth, 

Grid  Spacing 

Index 

meters 

Az.  meters 

1 

0 

2 

62 

61 

63 

62 

0.5 

186 

123.5 

187 

124 

1 

197 

134 

198 

135 

2 

208 

155 

209 

157 

4 

234 

257 

The  grid  spacing  in  the  range  direction  is  the  same  as 
in  Table  1  except  that  the  region  where  Ar  =  20  m  is 
extended  to  31*6  km. 
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Fig*  15*  Top:  Observed  25  Hz  airgun  signal  level  versus  range  found 
using  the  6.4  second  data  window  shown  in  Figure  12.  The  plot  is 

corrected  for  cylindrical  spreading. 
Middle:  25  Hz  parabolic  equation  output  at  9  m  water  depth  for  the 

acoustic  model  of  Figure  13. 

Bottom:  Parabolic  equation  output  for  same  model  as  used  above  with  the 

addition  of  a  200  m  thick  sediment  trough  situated  between  7  and  22  km 
range. 
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It  was  believed  that  the  airgun  operated  at  9  m  water  depth.  Using 
the  model  of  Figure  13,  the  acoustic  wave  field  at  9  m  depth  was 
calculated  for  25  Hz  (Figure  15,  middle).  The  result  of  the  calculation 
showed  fluctuating  signal  levels  near  the  source  which  indicated  the 
presence  of  several  propagating  modes.  The  smoothness  of  the  curve 
beyond  a  range  of  10  km  suggested  that  only  one  mode  existed  in  this 
region.  The  calculated  propagation  loss  between  15.4  and  29.5  km  was 
0.20  dB/km. 

In  both  the  observed  data  and  the  computed  signal  levels  there  was 
a  transition  frdm  multimode  to  single  mode  propagation.  However,  the 
range  where  this  transition  occurred  was  different  in  both  cases.  This 
range  mismatch  was  probably  due  to  the  use  of  an  attenuating  sub-bottom, 
a  point  which  is  discussed  later. 

In  an  effort  to  better  fit  the  measured  propagation  loss  rate,  a 
200  m  thick  sediment  trough  was  added  to  the  model  between  7  and  22  km. 
A  sediment  trough  could  represent  a  model  for  a  sediment  filled 
Pleistocene  valley  or  stream  bed  [King  and  MacLean,  1970].  The  computed 
propagation  loss  for  this  model  is  shown  in  Figure  15  (bottom).  The 
trough  extended  the  range  of  multimode  propagation  as  shown  by  the 
wiggliness  at  the  left  end  of  the  trough,  but  did  not  provide  a  better 
fit  to  the  data. 

A  sediment  velocity  gradient  of  1  s  added  to  the  above  sediment 
trough  model  produced  an  insignificant  change  in  the  propagation  loss 
results.  The  finding  that  a  velocity  gradient  in  the  sediments  did  not 
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influence  the  propagation  loss  indicated  that  most  of  the  energy  at  25 
Hz  was  trapped  in  the  water  column. 


Comparison  of  data  to  model  at  10  Hz 

The  10  Hz  spectrum  level  versus  range  is  plotted  in  Figure  16 
(top).  Values  for  ranges  between  15.4  and  27.5  km  are  displayed  and  the 
propagation  loss  was  0.14  dB/km  compared  to  the  -0.17  dB/km  loss  found 
at  25  Hz.  At  10  Hz  there  was  less  variation  in  signal  level  from  shot 
to  shot  than  at  25  Hz. 

Figure  16  also  shows  the  calculated  wavefield  at  10  Hz  and  at  a 
depth  of  9  m  in  the  water  column  plotted  versus  range.  The  calculated 
propagation  loss  between  15.4  and  27.5  km  of  -0.43  dB/km  does  not  fit 
the  data.  The  possibility  that  this  propagation  loss  rate  is  controlled 
by  the  sediment  column  is  suggested  by  the  observation  that  at  10  Hz, 
signals  propagating  at  1460  m/ s  have  a  wavelength  of  146  m  which  is 
twice  the  average  water  depth  of  about  70  m  and  greater  than  the 
combined  water  and  sediment  thickness  of  91  m.  In  this  case,  the  first  / 
mode  was  slightly  below  the  cutoff  frequency  and  the  calculated  signal 
level  rapidly  decayed  within  the  first  10  km  (Figure  16).  Note  that  for 
a  70  m  water  layer  with  velocity  1452  m/ s  over  a  sediment  half-space 
with  velocity  1600  m/s,  the  cutoff  frequency  is  12.3  Hz.  The  presence 
of  the  2000  m/ s  acoustic  basement  lowers  this  cutoff  frequency  slightly. 

The  model  sediment  layer  was  thin  in  comparison  with  the  acoustic 


PROP.  LOSS  y, 
+0 . 14  dB/km  I"  o 


59 


10  Hz 


| ! - 1 - - * 

5  |  0  20  30 

RANGE  km 


Fig.  16.  From  top  to  bottom: 

10  Hz  signal  levels  found  using  a  6.4  second  window  and  corrected  for 

cylindrical  spreading; 

Parabolic  equation  output  for  the  geoacoustic  model  shown  in  Figure  13 

at  10  Hz; 

10  Hz  output  for  above  model  with  200  m  thick  sediment  trough  from  7  to 

22  km  range; 

10  Hz  output  for  sediment  trough  model  with  Is1  velocity  gradient. 
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wavelength,  and  possible  shear  conversions  at  the  sediment-basement 
interface  [Vidmar,  1980]  which  cause  additional  propagation  loss  were 
neglected  by  the  parabolic  equation  used  here. 

Adding  the  200  m  thick  sediment  trough  between  7  and  22  km  gave  a 

propagation  loss  of  -0.15  dB/km  (Figure  16).  Placing  a  1  s  ^  gradient 
in  the  sediment  trough  model  resulted  in  a  realistic  propagation  loss  of 
0.12  dB/km  (Figure  16,  bottom)  but  the  variance  in  the  model  signal 
level  was  much  greater  than  that  of  the  data.  It  was  obvious  that,  as  a 
consequence  of  the  shallowness  of  the  water  with  respect  to  the  10  Hz 
wavelength,  the  sediments  controlled  the  propagation  loss  rate. 
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V.  DISCUSSION  ON  SOURCES  OF  MODELING  ERROR 

In  the  acoustic  modeling  presented  here,  several  approximations 
were  made  which  deserve  attention  as  likely  sources  of  error*  The 
parabolic  equation  is  limited  to  treating  a  narrow  band  of  propagating 
rays.  The  Epstein  layer  is  an  approximation  to  the  acoustic  interface 
and  causes  a  phase  error  in  the  reflected  wavefield.  Applying  the 
Gaussian  source  pulse  at  an  interface  may  cause  errors.  The  use  of  an 
attenuating  sub-bottom  is  not  physically  real  and  may  affect  propagation 
loss.  These  four  factors  are  discussed  below. 

Earlier  it  was  mentioned  that  equation  [7]  is  known  as  the  15°  or 
narrow  bandwidth  approximation.  In  the  model  for  propagation  on  the 
Nova  Scotian  continental  slope,  the  critical  angle  at  the  water- sediment 
interface  was  65.2°  .  Coincidently ,  the  same  angle  was  obtained  for 
critical  incidence  at  the  water-sediment  interface  in  the  Scotian  shelf 
model.  Since  most  of  the  propagating  modes  of  the  computed  wave  fields 
were  associated  with  postcritical  bottom  reflections,  the  incidence 
angles  associated  with  those  modes  ranged  from  65°  to  90°.  This  spread 
of  25°  was  thus  greater  then  the  15°  limit  of  equation  [7].  In  general, 
the  error  arising  from  using  the  15°  approximation  to  propagate  a  wider 
range  of  incident  rays  or  corresponding  phase  velocities  causes  a  shift 
in  the  positions  of  convergence  zones  [Brock  et  al.,  1977]. 

The  phases  of  the  Epstein  layer  and  Rayleigh  reflection 
coefficients  were  quite  different,  as  is  shown  in  Figures  6  and  14. 
This  phase  difference  was  unavoidable  since  the  Epstein  layer  "reflects” 
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rays  by  forcing  them  to  turn.  Since  the  wave  field  computed  from 
equation  [7]  is  a  standing  wave  interference  pattern,  the  incorrect 
phase  caused  by  using  Epstein  layers  can  be  expected  to  cause 
significant  shifts  in  the  position  of  local  field  values.  Choosing  the 
Epstein  layer  parameters  so  that  the  magnitude  of  the  Epstein  layer 
reflection  coefficient  matched  that  of  the  Rayleigh  interface,  however, 
ensured  that  the  energy  trapped  in  the  water  layer  was  identical  to  that 
when  acoustic  interfaces  are  included  in  the  modeling.  Thus,  one  would 
expect  the  propagation  loss  rates  found  over  long  distances  to  be 
similar  for  models  employing  Epstein  layers  and  models  using  true 
interfaces. 

There  is  also  some  inaccuracy  introduced  by  placing  the  Gaussian 
pulse  on  the  water-sediment  interface.  The  Gaussian  source  pulse  [23] 
was  selected  from  a  variety  of  possible  source  functions  because  it  is 
simple  to  compute.  An  alternative  would  be  to  use  as  initial  data  the 
normal  mode  wavefield  at  a  short  range  from  the  source  [Wood  and 
Papadakis,  1980].  It  seems  probable  that  the  symmetric  source  function 
[23]  did  not  accurately  model  the  source  on  the  interface.  The  water- 
sediment  interface  itself  was  modeled  as  a  smoothly  varying  transition 
zone.  Placing  a  Gaussian  pulse  on  such  a  transition  region  should  at 
least  be  more  physically  reasonable  than  locating  the  pulse  on  a 
discontinuity. 

The  use  of  an  attenuating  sub-bottom  causes  higher  order  modes 
associated  with  rays  at  precritical  bottom  incidence  to  be  prematurely 
damped  out.  This  result  may  be  seen  in  Figure  15  by  comparing  the 
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observed  and  predicted  25  Hz  signal  levels.  The  higher  order  modes  in 
the  observed  signal  level  are  present  out  to  about  25  km.  In  the 
predicted  signal  levels,  the  higher  order  modes  are  not  visible  beyond 
12  km.  In  this  example,  improper  modeling  of  the  bathymetry  and 
sediment  properties  may  also  be  responsible  for  the  premature  decay  of 
the  higher  modes.  However,  it  is  likely  that,  at  close  ranges,  use  of 
an  attenuating  sub-bottom  causes  erroneous  propagation  loss  when 
precritical  bottom  reflections  are  significant. 


64 


VI.  SUMMARY 


An  implicit  finite  difference  algorithm  for  the  parabolic  equation 
based  on  a  variable  grid  spacing  in  the  depth  as  well  as  range  direction 
was  developed.  The  algorithm,  when  implemented  with  a  variable  grid 
spacing,  gave  almost  the  same  results  as  a  constant  fine  grid  spacing  in 
one  third  as  much  computer  time. 

Acoustic  interfaces,  characterized  by  a  velocity  and  a  density 
contrast,  were  modeled  by  using  Epstein  layers.  The  Epstein  layer 
parameters  were  selected  by  fitting  the  amplitude  of  the  theoretical 
Epstein  layer  reflectivity  [ Brekhovshkikh ,  1980]  to  the  acoustic 
Rayleigh  reflectivity.  The  use  of  Epstein  layers  provided  an  interface¬ 
like  aspect  to  the  acoustic  models.  It  was  necessary,  however,  to 
specify  Epstein  layer  thicknesses  roughly  proportional  to  the  acoustic 
wavelength,  making  the  Epstein  layer  difficult  to  use  for  low  frequency 
modeling  of  interfaces  separated  by  a  thin  layer. 

Propagation  loss  curves  were  computed  using  this  algorithm  for  shot 
lines  on  the  Scotian  slope  and  shelf.  The  modeling  results  were 
compared  with  data  recorded  by  hydrophones  positioned  on  the  ocean 
bottom.  On  the  Scotian  slope,  the  observed  -0.19  dB/km  propagation  loss 
at  25  Hz  could  not  be  replicated.  The  observed  negative  propagation 
loss  may  be  partly  caused  by  compression  of  the  signals  by  the  recording 


instrumentation. 
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On  the  Scotian  shelf,  in  water  about  70  m  deep,  the  modeling  could 
not  explain  the  observed  -0.17  dB/km  propagation  loss  at  25  Hz.  At  this 
frequency,  the  sub-bottom  acoustic  parameters  did  not  greatly  effect  the 
modeled  propagation  loss.  This  result  was  in  agreement  with  the 
observation  that  25  Hz  was  well  above  the  cutoff  frequency  for  normal 
mode  propagation  in  the  water  column.  On  the  other  hand,  10  Hz  was 
close  to  the  cutoff  frequency  and  the  observed  0.14  dB/km  propagation 
loss  at  10  Hz  was  simulated  by  adjusting  the  properties  of  the 
sediments.  Seismic  reflection  profiling  on  the  slope  and  shelf  would  be 
helpful  for  determining  the  depth  to  acoustic  basement  and  extent  of 
lateral  sediment  changes  needed  to  improve  the  geoacoustic  models. 
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APPENDIX  A.  EPSTEIN  LAYER  REFLECTIVITY  FUNCTION 


The  parameters  of  equation  [31]  are  here  described  in  detail. 
S  =  2k^/m  is  the  relative  thickness  of  the  Epstein  layer.  6  is  the 
angle  of  incidence  of  the  plane  wave  upon  the  Epstein  layer.  M  is  the 
constant  of  equation  [28]  which  is  related  to  the  density  contrast 
across  the  interface  for  which  the  Epstein  layer  is  an  approximation. 
N  is  the  constant  that  is  determined  by  the  velocity  contrast  across  the 
interface.  a  ,  8,  and  Y  of  equation  [31]  are  specified  by  the  following 
equations: 

A  =  Ret  1(1  -  4S2M)1/2] 

2 

B  =  Im[  1(1  -  4S2M)1/2] 

2 

1  2  1/2 

a  =_  +  A  +  (iS/2)[cos8  -  (cos  8  -  N)  ]  +  iB 

2 

8  =  i  +  A  +  (iS/2)[cose  +  (cos20  -  N)^2]  +  iB 
2 

Y  =  1  +  iS  cos  9 

Brekhovskikn  [1980]  derived  the  above  expressions  by  relating  the 
solution  of  the  reduced  wave  equation  [2]  to  the  solution  of  the 
hypergeometric  equation. 
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The  Fortran  program  given  on  the  next  page  computes  ^(z+1)  for 
Re ( z)  >  -11/2.  For  values  of  z  outside  this  range,  continuation 
formulas  may  be  used.  Here  a  seven  term  series  expansion  due  to  Lanczos 
is  implemented  (Luke,  1980).  A  typical  error  for  In  T(z)  is  on  the 

order  of  10-7.  If  z  =  x  +  iy,  then  eZ  =  eX(cos(y)  +  isin(y)).  This 
relation  may  be  used  if  complex  exponentiation  is  not  provided  on  the 
user's  computer. 
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FUNCTION  CGAM(CZ) 

C  CGAM  IS  THE  GAMMA.  FUNCTION  FOR  COMPLEX  ARGUMENT  Z 
C  FOR  INPUT  COMPLEX  Z,  OUTPUT  IS  GAMMA (Z  +1). 

C 

C  REQUIRE  REAL ( Z )  >  -11/2 

C  HOWEVER,  IT  IS  STRONGLY  RECOMMENDED  THAT  YOU  USE  CONTINUATION 
C  FORMULAS  FOR  REAL(Z)  <  0  BECAUSE  OF  THE  SINGULARITES  IN 
C  GAMMA (Z+l)  WHEN  Z  =  -1 ,  -2,  -3,  ETC. 

C 

C  THIS  IS  THE  LANCZOS  EXPANSION  FROM 

C  Y.L.  LUKE,  THE  SPECIAL  FUNCTIONS  AND  THEIR  APPLICATIONS, 

C  ACADEMIC,  N.Y.,  1969,  VOL.  1  &  2. 

IMPLICIT  COMPLEX(C) 

DIMENSION  G(0: 6) 

C  THERE  ARE  (NORD  +1)  G  CONSTANTS 
5  NORD  =  6 

G( 0)  =  41.624436916439068 

G( 1)  =  -51.224241022374774 
G( 2)  =  11.338755813488977 

G(3)  =  -0.747732687772388 

G( 4)  =  0.008782877493061 

G( 5)  =  -0.000001899030264 

G(6)  =  .  0.000000001946335 

•PI  =  3.141592654 
E  =  2.718281828 
PI2SQ  =  SQRT(2 .0  *  PI) 

Cll  =  CZ  +  11.0  /  2.0 

C2  =  CZ  +  0.5 

C2  =  Cll  **  C2  *  PI2SQ 

C2  =  C2  *  E  **  (Cll  *  (-1.0)) 

CH  =  (1.0, 0.0) 

10  CSUM  =  (0.0, 0.0) 

FOR  1=0,  NORD 
.  CSUM  =  CSUM  +  CH  *  G(I) 

.  RIX  =  I 
.  RIX1  =1+1 

20  .  CH  =  CH  *  (CZ  -  RIX)  /  (CZ  +  RIX1) 

END  FOR 

CGAM  =  C2  *  CSUM 

RETURN 

END 
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APPENDIX  B.  RAYTRACING  ON  THE  CONTINENTAL  SLOPE 

This  appendix  describes  raytrace  modeling  of  acoustic  signals 
recorded  on  the  continental  slope  off  Nova  Scotia*  The  modeled  signals 
are  the  same  as  those  modeled  using  the  parabolic  equation  as  described 
in  Chapter  2.  The  SUS  charge  line  was  downslope  from  the  telemetering 
ocean  bottom  seismometer  (TOBS  3)  on  an  average  slope  of  3*5°.  At  TOBS 
3  the  water  depth  was  1301  m,  while  at  the  furthest  shot  the  water  depth 
was  over  2800  m.  The  shots  were  1.1  oz  (0.031  kg)  SUS  charges  detonated 
at  a  depth  of  18  m.  Figure  2  shows  the  location  of  the  SUS  charge  line 
on  the  continental  slope,  as  well  as  the  position  of  two  sonic  velocity 
versus  depth  profiles  of  the  water  column  (Figure  3). 

For  this  snot  line  the  water  was  sufficiently  deep  to  allow 
individual  bottom  bounce  multiple  arrivals  of  the  water  wave  to  be 
readily  identified.  Yet,  the  ocean  was  shallow  enough  so  that  beyond  a 
few  km  range,  the  most  energetic  arrivals  were  all  reflections  off  the 
bottom.  In  addition,  the  average  3.5°  slope  (local  slopes  were  as  much 
as  18°)  caused  the  angle  of  incidence  of  a  downslope-propagating  ray  to 
increase  upon  each  bottom  bounce;  the  incidence  angle  of  an  upslope- 
propagating  ray  would  have  decreased. 

Ray  tracing  is  a  useful  tool  for  this  problem  because  it  enables 
the  observer  to  visualize  actual  wave  propagation  effects.  Water  wave 
arrivals  recorded  by  TOBS  3  and  the  ray  traces  for  one,  two,  and  three 
bottom  multiple  arrivals  are  pictured  in  Figures  17  to  19.  In  these  ray 
trace  figures,  the  rays  emanate  from  the  position  of  TOBS  3  in  0.5 
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degree  increment s •  A  composite  of  sonic  velocity  stations  4  and  5  is 
used  for  the  water  velocity  versus  depth  function.  The  bottom  reflects 
all  incident  rays  and  is  slightly  smoothed  to  eliminate  errant  rays 
caused  by  irregularities  in  the  bathymetry  profile. 

For  single  bottom  bounce  water  wave  arrivals  (Figure  17),  the 
shape  ot  the  ocean  bottom  is  the  most  significant  factor  affecting 
propagation  loss.  A  convex  bottom  near  9  km  creates  a  small  shadow  zone 
at  15  km  range  while  a  convex  bottom  near  14  km  results  in  a  shadow  zone 
from  24  to  30  km  (up  to  the  limit  of  the  plot),  A  concave  bottom  near 
13  km  range  produces  a  region  of  high  intensity  between  17  and  18  km 
which  shows  up  in  the  data  between  16  and  17  km  range.  This  discrepancy 
is  probably  caused  by  slight  errors  in  modeling  the  bathymetry  and  water 
column  velocity.  Another  possible  source  of  error  is  neglect  of 
azimuthal  dependence  in  the  raytracing. 

The  two  and  three  bottom  bounce  multiple  arrivals  indicated  in 
Figures  18  and  19  largely  consist  of  postcritical  reflections  since 
suffering  more  than  one  precritical  reflection  greatly  reduces  the 
amplitude  of  the  arrival.  Beyond  13  km  surface  range,  the  increased 
amplitude  of  the  second  multiple  (Figure  18)  is  due  to  the  second  bounce 
reaching  or  exceeding  critical  incidence  of  about  65°  (incident  angles 
are  measured  from  the  normal  to  the  bottom),  starting  from  around  9.5  km 
range  on  the  bottom.  For  this  discussion,  bottom  bounces  are  counted 
from  the  first  bounce  downsiope  of  the  ray  source  on  the  bottom. 
Between  13  and  18  km  range,  the  second  bounce  arrival  is  additionally 
intensified  by  bottom  focusing. 
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Careful  examination  reveals  some  increase  in  the  amplitude  of  the 
third  multiple  arrivals  (Figure  19)  between  14  and  16  km,  due  to  bottom 
focusing.  However,  only  past  22  km  do  the  third  multiples  become 
significant,  a  consequence  of  their  becoming  critically  incident  at 
these  ranges.  The  family  of  post  critical  third  bounce  arrivals  after 
22  km  originates  from  the  second  multiples  seen  between  13  and  18  km. 


Fig.  17.  One  bottom  bounce  ray  paths  from  TOBS  3.  The  shape  of  the 
ocean  floor  creates  a  shadow  zone  beyond  24  km.  The  solid  line  precedes 
the  single  bottom  bounce  arrivals  * 
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Fig,  180  Double  bottom  bounce  ray  paths  from  TOBS  3.  The  second  bounce 
achieves  critical  incidence  of  65°  at  around  9*5  km  range  on  the  bottom. 
The  solid  line  precedes  the  second  bottom  bounce  arrivals. 
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APPENDIX  C.  HEAD  WAVE 

A  phase  reversal  at  short  range  was  observed  in  the  onset  of 
explosive  signals  received  by  TOBS  3  from  large  1.8  lb  (0.82  kg)  SUS 
charges.  The  large  SUS  charge  line  was  along  the  small  SUS  charge  line 
of  Figure  2  and  extended  upslope  from  the  receiver.  The  shots  were  set 
off  at  18  m  depth  and  were  recorded  in  water  1301  m  deep. 

Inspection  of  the  first  breaks  in  the  first  arrivals  reveals  a 
phase  reversal  (Figure  20)  at  shots  with  ranges  of  2.0  and  3.8  km 
upslope  and  downslope  from  TOBS  3,  respectively.  It  is  probable  that 
the  phase  change  indicates  the  onset  of  a  critically  incident  head  wave. 
The  developments  below  assume  that  the  upslope  and  downslope  phase 
reversals  mark  rays  that  impinge  upon  the  bottom  at  the  same  incidence 
angle  (Figure  21).  This  incidence  angle  is  the  critical  angle,  from 
which  the  sediment  velocity  may  be  computed  if  the  water  velocity  is 
known. 

By  assuming  a  direct  water  borne  path  from  a  shot  near  the  surface 
to  the  receiver  on  the  ocean  bottom,  the  ray  parameter  corresponding  to 
the  travel  time  of  a  first  arrival  may  be  calculated.  The  water 
velocity  structure  on  the  left  side  of  Figure  3,  extrapolated  linearly 
down  to  1301  m,  is  used  for  this  calculation.  For  the  water  velocity 
structure  and  the  shot  ranges  under  consideration,  it  happens  that  each 
travel  time  corresponds  to  a  unique  ray  parameter.  The  ray  parameter  is 
computed  from  the  following  arrangement  of  a  well  known  formula ; 
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Fig.  20.  Above:  Large  SUS  charge  shots  originating  upslope  from  TOBS  3 
Below:  Shots  originating  downslope  from  TOBS  3. 

The  triangles  denote  the  phase  reversal  in  the  first  break. 
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Fig.  21.  Top:  Geometry  for  ray  paths  occurring  at  65.2°  critical 
°  incidence. 

Bottom:  Bottom  angles  a  and  3  are  computed  from  their  corresponding  ray 

parameters  and  the  bottom  water  velocity,  1488  m/s.  <t>  is  the  bottom 
slope  in  the  vicinity  of  TOBS  3  and  may  be  computed  from  a  and  8 . 

Using  <t>  and  a  or  8  ,  the  65.2°  incident  angle  with  respect  to  the  bottom 
is  determined. 
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where  z^  *  shot  depth,  D  *  bottom  depth,  T  =  travel  time,  n  =  l/c(z), 
and  p  **  p(T)  =  ray  parameter.  Equation  [33]  may  be  solved  for  p  by 
providing  an  initial  guess  for  p  and  using  the  secant  method  or  other 
root  finding  technique.  The  integrand  is  numerically  integrated. 

Once  the  ray  parameter  is  found  for  the  rays  corresponding  to  the 
upslope  and  the  downslope  phase  reversal,  the  ocean  bottom  slope  and 
bottom  incident  angle  associated  with  the  phase  reversal  may  be 
computed,  as  shown  in  Figure  21.  This  bottom  incident  angle  is  the 
critical  angle. 

The  nature  of  the  phase  reversal  is  investigated  further  by 
examining  the  first  break  amplitudes.  The  ray  parameter  for  each  shot 
is  obtained  using  equation  [33].  From  each  ray  parameter,  the  bottom 
incident  angle  and  horizontal  range  are  calculated.  First  break 
amplitudes  for  shots  downslope  from  the  receiver  versus  bottom  incident 
angle  and  horizontal  range  are  plotted  in  Figure  22.  The  expected 
amplitude  decay  due  to  spherical  spreading  is  also  indicated  in  the 
figure.  The  failure  of  shots  closer  than  1600  m  to  obey  spherical 
spreading  suggests  that  the  near  shot  first  breaks  saturated  the 
recording  amplifier.  The  phase  reversal  occurs  abruptly  just  past  65° 
bottom  incidence,  signaling  the  onset  of  the  head  wave.  Beyond  65°,  the 
first  arrival  is  no  longer  the  direct  wave  in  the  water  so  that  the 
computed  incident  angles  and  ranges  no  longer  apply.  This  fact  is 
illustrated  in  the  lower  part  of  Figure  22.  Shown  is  the  difference 
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between  horizontal  ranges  based  on  equation  [33]  ray  parameters  and 
horizontal  ranges  obtained  from  shipboard  navigation.  The  sudden  drop 
off  beyond  65°  incidence  indicates  that  the  computed  ray  parameters  no 
longer  describe  the  first  arrival,  which  is  no  longer  a  water  borne 
direct  wave.  The  same  results  are  obtained  by  looking  at  the  first 
break  amplitudes  for  shots  upslope  of  the  TOBS  (Figure  23).  In  this 
case,  the  head  wave  onset  is  indicated  by  a  sharp  upward  swing  of  the 


range  difference  plot. 
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Fig.  22.  Above:  First  break  amplitudes  versus  computed  bottom  incident 
angle  and  horizontal  range  of  raypath.  The  shots  are  downslope  from  the 
receiver.  Solid  line  is  loss  caused  by  spherical  spreading,  normalized 

on  the  last  shot  before  the  phase  reversal. 
Below:  Difference  between  the  computed  horizontal  range  and  the  range 

from  shipboard  navigation  versus  calculated  bottom  incident  angle. 
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Fig.  23.  Above:  First  break  amplitudes  versus  computed  bottom  incident 
angle  and  horizontal  range  of  raypath  for  shots  upslope  from  the 
receiver.  Solid  line  is  spherical  spreading  loss,  normalized  on  the 

last  shot  before  the  phase  reversal. 
Below:  Difference  between  the  computed  horizontal  range  and  the  range 

from  shipboard  navigation  versus  calculated  bottom  incident  angle. 
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APPENDIX  D.  SEDIMENT  VELOCITY  FROM  APPARENT  VELOCITIES 


The  apparent  velocities  measured  from  the  leading  and  trailing 

edges  of  the  water  wave,  as  shown  in  Figure  12,  may  be  used  to  estimate 

the  ocean  bottom  sediment  velocity.  The  method  for  doing  so  was 

developed  by  Sutton  and  Maynard  [  1 97 1 J  and  described  by  Houtz  [  1980]. 

The  technique  is  recapitulated  here. 

Consider  a  single  shot.  Let  v  be  the  water  velocity  and  v  be  the 

w  s 

bottom  velocity.  Let  t^  be  the  travel  time  of  the  direct  wave  which  is 

at  the  leading  edge  of  the  water  wave.  The  travel  time  of  the  back  end 

of  the  water  wave  is  t  .  This  arrival  corresponds  to  the  ray  traveling 

c 

at  critical  incidence  9  with  respect  to  the  bottom  (Figure  24).  It  is 

c 

assumed  that  the  bottom  sediment  velocity  is  greater  than  the  water 
velocity.  The  apparent  velocity  of  the  terminal  arrival  is  vg  *  v^t^/t^. 

From  Figure  24  the  length  AB  =  v^t^,  and  AB  =  v^t^sin  9  ^ .  Using  the 
above  relations,  one  may  write: 


sin  9 

c 

From  Snell's  law, 


v  = 
s 


w 


sin  9 


Replacing  sin  0^  yields  the  result: 


v  =  v  /  v 
s  w  g 


[34] 


[35] 


[36] 


\ 
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From  Figure  12  v  =  1.42  km/s  and  v  =  1.32  km/s.  Using  equation  [36], 
w  § 

V  «  1.53  km/s.  This  velocity  is  a  bit  lower  than  Scotian  shelf  values 
s 

of  1.57  to  1.66  km/s  measured  by  McKay  and  McKay  [1982J.  The  slope 
correction  given  by  Houtz  [1980J  is  negligible  for  this  example.  One 
source  of  error  may  be  the  assumption  of  constant  water  velocity.  The 
near  surface  gradient  in  the  water  velocity  versus  depth  shown  in  Figure 
13  hastens  the  arrival  of  the  critically  reflected  wave  so  that  sin 
of  equation  [34]  is  inflated.  The  increased  sin  produces  a  lower 
calculated  sediment  velocity  v  ,  as  shown  in  equation  [35]. 


Fig.  24.  Raypath  at  critical  incidence  0  .  The  last  arrival  in  the 
water  wave  travels  at  critical  incidence#  The  total  distance  traveled 
is  the  hypotenuse  of  the  large  triangle,  or  v  t  (see  text). 
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